Measurement exploratory data analysis¶

In [ ]:
DEVICES = [
    'beaglebone-fan',
    'beaglebone-compressor',
    'beaglebone-pump',
    'beaglebone-refrigerator',
    'mafaulda-a',
    'mafaulda-b'
]

DEVICE = DEVICES[0]
T_WAVEFORM = 5  # (1 = MaufaulDa, 5 = others)
T_SEC = T_WAVEFORM
NFFT = 2**10   # (2**10 = MaufaulDa, 2**14 = others)
F_LIMIT = None # (3000 = MaufaulDa, None = others)
In [ ]:
import os
import pandas as pd
import numpy as np
import matplotlib.pylab as plt

from tabulate import tabulate
from IPython.display import Markdown, HTML
from tqdm.notebook import tqdm

from typing import List, Tuple
from scipy.signal import find_peaks, butter, lfilter
from tsfel.feature_extraction.features import fundamental_frequency

from zipfile import ZipFile
import sys
sys.path.append('../')
from vibrodiagnostics import mafaulda, selection, discovery, models


def beaglebone_measurement(filename: str, fs: int) -> Tuple[str, pd.DataFrame]:
    g = 9.81
    milivolts = 1800
    resolution = 2**12
    columns = ['x', 'y', 'z']
    ts = pd.read_csv(filename, delimiter='\t', index_col=False, header=None, names=columns)
        
    # Calculate amplitude in m/s^2 Beaglebone Black ADC and ADXL335 resolution (VIN 1.8V, 12bits)
    for dim in columns:
        ts[dim] = ts[dim] * (milivolts / resolution)  # ADC to mV
        ts[dim] = (ts[dim] / 180) * g                 # mV to m/s^2 (180 mV/g)
        ts[dim] -= ts[dim].mean()

    ts['t'] = ts.index * (1 / fs)
    ts.set_index('t', inplace=True)
    return (os.path.basename(filename), ts, fs, ts.columns)  # last is feature columns


def beaglebone_dataset(filenames: List[str], fs: int) -> List[Tuple[str, pd.DataFrame]]:
    dataset = []
    for filename in filenames:
        name, ts, fs, cols = beaglebone_measurement(filename, fs)
        dataset.append((name, ts))
    return dataset


def lowpass_filter(data, cutpoint, fs, order=5):
    b, a = butter(order, cutpoint, fs=fs, btype='lowpass')
    y = lfilter(b, a, data)
    return y


def mafaulda_dataset(
        place=['ax', 'ay', 'az'],
        features_path =  '../../datasets/features_data/',
        mafaulda_path='../../datasets/MAFAULDA.zip',
        rpm=2500,
        lowpass_hz=10000):

    metadata_filename = os.path.join(features_path, selection.MAFAULDA_METADATA)
    faults = {
        'normal': 'normal',
        'imbalance': 'imbalance',
        'horizontal-misalignment': 'misalignment',
        'vertical-misalignment': 'misalignment',
        'overhang-cage_fault': 'cage fault',
        'underhang-cage_fault': 'cage fault',
        'underhang-ball_fault': 'ball fault',
        'overhang-ball_fault': 'ball fault',
        'overhang-outer_race': 'outer race fault',
        'underhang-outer_race': 'outer race fault',
    }

    metadata = pd.read_csv(metadata_filename, index_col='filename')
    metadata.reset_index(inplace=True)
    metadata = models.fault_labeling(metadata, faults)
    files = pd.DataFrame()
    # Worst severity and mid rpm
    for name, group in metadata[metadata['rpm'] >= rpm].groupby(by='fault', observed=False):
        files = pd.concat([
            files,
            group[
                group['severity_level'] == group['severity_level'].max()
            ].sort_values(by='rpm', ascending=True).head(1)
        ])
    ordering = {
        'normal': 0,
        'misalignment': 1,
        'imbalance': 2,
        'cage fault': 3,
        'ball fault': 4,
        'outer race fault': 5,
    }
    source = ZipFile(mafaulda_path)
    dataset = len(files) * [0]
    for index, file in files.iterrows():
        ts = mafaulda.csv_import(source, file['filename'])
        ts = ts[place]
        ts.columns = ts.columns.str.extract(r'(\w)$')[0]
        for axis in ts.columns:
            ts[axis] = lowpass_filter(ts[axis], lowpass_hz, file['fs'])
        pos = ordering[file['fault']]
        dataset[pos] = ((file['fault'] + ' (' + file['filename'] +')', ts))

    return dataset

Load dataset

In [ ]:
if DEVICE == 'beaglebone-fan':
    Fs = 2500
    path = '../../inspections/fan/'
    files = [
        '1_still.tsv', '2_still.tsv', '3_still.tsv',
        '1_up.tsv', '2_up.tsv', '3_up.tsv',
        '1_down.tsv', '2_down.tsv', '3_down.tsv'
    ]
    files = [os.path.join(path, name) for name in files]
    DATASET = beaglebone_dataset(files, Fs)

elif DEVICE == 'beaglebone-compressor':
    Fs = 2500
    path = '../../inspections/datacentres/shc3/'
    files = [
        'k3_1.tsv', 'k3_2.tsv', 'k3_3.tsv', 'k3_4.tsv',
        'k5_1.tsv', 'k5_2.tsv', 'k5_3.tsv', 'k5_4.tsv'
    ]
    files = [os.path.join(path, name) for name in files]
    DATASET = beaglebone_dataset(files, Fs)

elif DEVICE == 'beaglebone-pump':
    Fs = 2500
    path = '../../inspections/bvs/'
    files = [
        'bvs_1_hore.tsv', 'bvs_2_hore.tsv' 
        #, 'bvs_3_motor.tsv', 'bvs_4_motor.tsv'
    ]    
    files = [os.path.join(path, name) for name in files]
    DATASET = beaglebone_dataset(files, Fs)

elif DEVICE == 'beaglebone-refrigerator':
    Fs = 2500
    path = '../../inspections/home-refrigerator/'
    files = [
        'ch1.tsv', 'ch2.tsv', 'ch3.tsv', 'ch4.tsv', 'ch5.tsv'
    ]    
    files = [os.path.join(path, name) for name in files]
    DATASET = beaglebone_dataset(files, Fs)

elif DEVICE == 'mafaulda-a':
    Fs = mafaulda.FS_HZ
    DATASET = mafaulda_dataset(place=['ax', 'ay', 'az'])

elif DEVICE == 'mafaulda-b':
    Fs = mafaulda.FS_HZ
    DATASET = mafaulda_dataset(place=['bx', 'by', 'bz'])
In [ ]:
for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    ts.info()
    print()

1_still.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 38200 entries, 0.0 to 15.2796
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       38200 non-null  float64
 1   y       38200 non-null  float64
 2   z       38200 non-null  float64
dtypes: float64(3)
memory usage: 1.2 MB

2_still.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 38200 entries, 0.0 to 15.2796
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       38200 non-null  float64
 1   y       38200 non-null  float64
 2   z       38200 non-null  float64
dtypes: float64(3)
memory usage: 1.2 MB

3_still.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 38200 entries, 0.0 to 15.2796
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       38200 non-null  float64
 1   y       38200 non-null  float64
 2   z       38200 non-null  float64
dtypes: float64(3)
memory usage: 1.2 MB

1_up.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 38200 entries, 0.0 to 15.2796
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       38200 non-null  float64
 1   y       38200 non-null  float64
 2   z       38200 non-null  float64
dtypes: float64(3)
memory usage: 1.2 MB

2_up.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 38200 entries, 0.0 to 15.2796
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       38200 non-null  float64
 1   y       38200 non-null  float64
 2   z       38200 non-null  float64
dtypes: float64(3)
memory usage: 1.2 MB

3_up.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 38200 entries, 0.0 to 15.2796
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       38200 non-null  float64
 1   y       38200 non-null  float64
 2   z       38200 non-null  float64
dtypes: float64(3)
memory usage: 1.2 MB

1_down.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 38200 entries, 0.0 to 15.2796
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       38200 non-null  float64
 1   y       38200 non-null  float64
 2   z       38200 non-null  float64
dtypes: float64(3)
memory usage: 1.2 MB

2_down.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 38200 entries, 0.0 to 15.2796
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       38200 non-null  float64
 1   y       38200 non-null  float64
 2   z       38200 non-null  float64
dtypes: float64(3)
memory usage: 1.2 MB

3_down.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 38200 entries, 0.0 to 15.2796
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       38200 non-null  float64
 1   y       38200 non-null  float64
 2   z       38200 non-null  float64
dtypes: float64(3)
memory usage: 1.2 MB

In [ ]:
for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    display(tabulate(ts.describe(), headers='keys', tablefmt='html'))
    ts.boxplot(grid=True)
    plt.show()

1_still.tsv

x y z
count38200 38200 38200
mean 4.53297e-15 3.95579e-15 3.23185e-15
std 0.874418 1.2639 0.0920418
min -1.83871 -2.52106 -0.348954
25% -0.5933 -1.108 -0.061552
50% 0.00545463 -0.00628912 0.0102986
75% 0.580259 1.11937 0.058199
max 1.96937 2.43663 0.345601
No description has been provided for this image

2_still.tsv

x y z
count38200 38200 38200
mean 8.66044e-16 6.06082e-15 2.94447e-15
std 2.56509 2.48179 0.322518
min -4.5757 -3.73282 -0.952338
25% -2.44414 -2.63111 -0.281733
50% 0.429887 0.0752644 0.00566968
75% 2.27405 2.56608 0.269122
max 4.04637 3.66779 0.987628
No description has been provided for this image

3_still.tsv

x y z
count38200 38200 38200
mean -8.10354e-15 -3.27184e-16 1.05735e-14
std 2.15156 2.64354 0.396043
min -3.02308 -4.12174 -0.786905
25% -2.13692 -2.58893 -0.379752
50% -0.244853 -0.0502045 0.00345146
75% 1.95857 2.56635 0.362704
max 3.92248 4.21293 0.841708
No description has been provided for this image

1_up.tsv

x y z
count38200 38200 38200
mean -2.19487e-17 1.16068e-15 -4.63229e-15
std 0.850544 1.19633 0.160765
min -2.0824 -2.62078 -0.818255
25% -0.381937 -0.704763 -0.0757992
50% 0.00126648 -0.010207 -0.00394865
75% 0.40842 0.756199 0.0679019
max 2.01308 2.57641 0.930109
No description has been provided for this image

2_up.tsv

x y z
count38200 38200 38200
mean -3.61707e-15 -5.37892e-15 3.77964e-16
std 2.54025 2.40738 0.506515
min -5.14075 -4.29866 -1.69449
25% -1.97932 -2.23894 -0.305381
50% 0.00854621 0.0123751 0.0778218
75% 1.99641 2.16789 0.341274
max 5.58894 4.61081 1.44298
No description has been provided for this image

3_up.tsv

x y z
count38200 38200 38200
mean 3.0704e-15 -3.60256e-15 1.23954e-15
std 1.69868 2.35114 0.38339
min -3.6383 -4.07792 -0.81214
25% -1.53068 -2.18586 -0.333137
50% -0.0697164 0.0415091 0.00216618
75% 1.34335 2.10123 0.337469
max 3.92997 4.23279 1.12783
No description has been provided for this image

1_down.tsv

x y z
count38200 38200 38200
mean -1.9047e-16 -9.52239e-15 -2.13758e-15
std 0.272355 0.355092 0.0934928
min -2.8127 -2.01469 -0.868043
25% -0.0344782 -0.0268261 -0.0297866
50% -0.0105281 -0.0028759 0.0181137
75% 0.0134221 0.0210743 0.0420639
max 3.3425 2.05684 0.664769
No description has been provided for this image

2_down.tsv

x y z
count38200 38200 38200
mean 4.42564e-15 4.22792e-16 2.48355e-15
std 0.910304 0.87873 0.242263
min -5.11725 -3.82467 -1.76085
25% -0.015856 -0.0405354 -0.0124892
50% 0.00809416 0.007365 0.0354112
75% 0.0320444 0.0313152 0.0833116
max 4.89393 3.9352 1.56822
No description has been provided for this image

3_down.tsv

x y z
count38200 38200 38200
mean -4.89177e-15 -9.63697e-16 -5.49685e-15
std 0.675041 0.924902 0.236818
min -3.52907 -3.8726 -3.24596
25% -0.0323384 -0.0405711 -0.0126823
50% 0.015562 0.00732926 0.0591683
75% 0.0395122 0.0552297 0.0831185
max 3.77574 4.22256 7.81903
No description has been provided for this image

Statistical tests

  • Normality test: Kolmogorov–Smirnov test
  • Normality visual test: Quantile-quantile plot on chosen recording
  • Stationarity test: Augmented Dickey–Fuller test
  • Stationarity visual test: Autocorrelation plot
In [ ]:
from statsmodels.tsa.stattools import adfuller
from statsmodels.api import qqplot
from scipy.stats import kstest

normality_tests = []
for name, ts in DATASET:
    for x in ts.columns:
        p_value = kstest(ts[x], 'norm').pvalue
        test = {'name': name, 'axis': x, 'p-value': p_value, 'not-normal': p_value < 0.05}
        normality_tests.append(test)

normality_tests = pd.DataFrame.from_records(normality_tests)
print(normality_tests.value_counts('not-normal'))
normality_tests.describe()
not-normal
True    27
Name: count, dtype: int64
Out[ ]:
p-value
count 2.700000e+01
mean 3.994415e-66
std 2.075559e-65
min 0.000000e+00
25% 0.000000e+00
50% 0.000000e+00
75% 0.000000e+00
max 1.078492e-64
In [ ]:
name, ts = DATASET[0]
fig, ax = plt.subplots(1, len(ts.columns), figsize=(10, 4))
for i, x in enumerate(ts.columns):
    qqplot(ts[x], line='45', ax=ax[i], marker='.', alpha=0.5)
    ax[i].set_title(f'Axis: {x}')

plt.tight_layout()
print(name)
plt.show()
1_still.tsv
No description has been provided for this image
In [ ]:
stationarity_tests = []
for name, ts in tqdm(DATASET):
    for x in ts.columns:
        result = adfuller(ts[x].loc[T_WAVEFORM:T_WAVEFORM+1])
        p_value = result[1]
        test = {
            'name': name,
            'axis': x,
            'statistic': result[0],
            'p-value': p_value,
            'stationary': p_value < 0.001
        }
        stationarity_tests.append(test)

stationarity_tests = pd.DataFrame.from_records(stationarity_tests)
print(stationarity_tests.value_counts('stationary'))
stationarity_tests['p-value'].describe()
  0%|          | 0/9 [00:00<?, ?it/s]
stationary
True     26
False     1
Name: count, dtype: int64
Out[ ]:
count    2.700000e+01
mean     1.942588e-04
std      9.026014e-04
min      0.000000e+00
25%      3.192522e-29
50%      1.201432e-20
75%      4.274780e-14
max      4.677476e-03
Name: p-value, dtype: float64
In [ ]:
name, ts = DATASET[0]
fig, ax = plt.subplots(1, len(ts.columns), figsize=(10, 4))
for i, x in enumerate(ts.columns):
    ax[i].acorr(ts[x], maxlags=50)
    ax[i].set_title(f'Axis: {x}')

plt.tight_layout()
print(name)
plt.show()
1_still.tsv
No description has been provided for this image

Time domain histogram

In [ ]:
for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    axis = ts.columns
    ax = ts[axis].hist(figsize=(15, 3), grid=True, bins=100, layout=(1, 3), edgecolor='black', linewidth=0.5)
    plt.show()

1_still.tsv

No description has been provided for this image

2_still.tsv

No description has been provided for this image

3_still.tsv

No description has been provided for this image

1_up.tsv

No description has been provided for this image

2_up.tsv

No description has been provided for this image

3_up.tsv

No description has been provided for this image

1_down.tsv

No description has been provided for this image

2_down.tsv

No description has been provided for this image

3_down.tsv

No description has been provided for this image

Time domain waveform

In [ ]:
for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    axis = ts.columns
    
    ax = ts[axis].plot(figsize=(20, 8), grid=True, subplots=True)
    for i, axname in enumerate(axis):
        ax[i].set_xlabel('Time [s]')
        ax[i].set_ylabel(f'Amplitude ({axname}) [m/s^2]')
    plt.show()               # plt.savefig('waveform.png')

1_still.tsv

No description has been provided for this image

2_still.tsv

No description has been provided for this image

3_still.tsv

No description has been provided for this image

1_up.tsv

No description has been provided for this image

2_up.tsv

No description has been provided for this image

3_up.tsv

No description has been provided for this image

1_down.tsv

No description has been provided for this image

2_down.tsv

No description has been provided for this image

3_down.tsv

No description has been provided for this image

Time domain waveform zoom detail

In [ ]:
for name, ts in DATASET:
    axis = ts.columns
    display(Markdown(f'**{name}**'))
    ax = (ts[axis].iloc[int(T_WAVEFORM*Fs):int(T_WAVEFORM*Fs)+Fs]
                  .plot(figsize=(20, 10), grid=True, subplots=True))
    
    for i, axname in enumerate(axis):
        ax[i].set_xlabel('Time [s]')
        ax[i].set_ylabel(f'Amplitude ({axname}) [m/s^2]')
        plt.show()      # plt.savefig('waveform_zoom.png')

1_still.tsv

No description has been provided for this image

2_still.tsv

No description has been provided for this image

3_still.tsv

No description has been provided for this image

1_up.tsv

No description has been provided for this image

2_up.tsv

No description has been provided for this image

3_up.tsv

No description has been provided for this image

1_down.tsv

No description has been provided for this image

2_down.tsv

No description has been provided for this image

3_down.tsv

No description has been provided for this image

Time domain waveform zoom - faults side by side

In [ ]:
fig, ax = plt.subplots(len(DATASET), 3, figsize=(12, 15), sharex=True)

for idx, df in enumerate(DATASET):
    name, ts = df
    columns = ts.columns
    ax[idx][1].set_title(name)
    ax[idx][0].set_ylabel('Amplitude [m/s^2]')

    for pos, axis in enumerate(columns):
        data = ts[axis].loc[T_WAVEFORM:T_WAVEFORM+0.3]
        ax[idx][pos].plot(data.index, data, linewidth=1, color='darkblue')
        ax[idx][pos].grid()
    

plt.tight_layout()
plt.show()
No description has been provided for this image
In [ ]:
def spectogram(x, debug=True):
    fig, ax = plt.subplots(figsize=(15, 4))
    cmap = plt.get_cmap('inferno')
    pxx, freqs, t, im = plt.specgram(
        x, NFFT=NFFT, Fs=Fs,
        detrend='mean',
        mode='magnitude', scale='dB',
        cmap=cmap, vmin=-60
    )
    fig.colorbar(im, aspect=20, pad=0.04)
    ax.set_xlabel('Time [s]')
    ax.set_ylabel('Frequency [Hz]')
    mafaulda.resolution_calc(Fs, NFFT)
    return freqs, pxx


def window_idx(t):
    return (Fs * t) // NFFT + 1


def spectrum_slice(freqs, Pxx, t):
    fig, ax = plt.subplots(2, 1, figsize=(20, 8))
    n = window_idx(t)

    dB = 20 * np.log10(Pxx.T[n] / 0.000001)
    ax[0].plot(freqs, dB)      # 1 dB = 1 um/s^2
    ax[0].grid(True)
    ax[0].set_xlabel('Frequency [Hz]')
    ax[0].set_ylabel('Amplitude [dB]')
    
    ax[1].plot(freqs, Pxx.T[n])
    ax[1].grid(True)
    ax[1].set_xlabel('Frequency [Hz]')
    ax[1].set_ylabel('Amplitude [m/s^2]')
    return n


def get_max_frequency(freqs, Pxx, i):
    max_freq = freqs[np.argmax(Pxx.T[i])]
    return max_freq


def get_peaks(freqs, Pxx, i, top=5):
    amplitudes = Pxx.T[i]
    peaks, _ = find_peaks(amplitudes, distance=3)

    fundamental = get_max_frequency(freqs, Pxx, i)
    f_top = freqs[peaks[np.argsort(amplitudes[peaks])]][::-top]
    y_top = np.sort(amplitudes[peaks])[::-top]

    return pd.DataFrame({
        'f': f_top,
        'y': y_top,
        '1x': f_top / fundamental 
    })


def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    b, a = butter(order, [lowcut, highcut], fs=fs, btype='band')
    y = lfilter(b, a, data)
    return y


def get_spectrograms(DATASET: List[pd.DataFrame], axis: str) -> list:
    spectrograms = []

    for name, ts in DATASET:
        base_freq = fundamental_frequency(ts[axis], Fs)
        display(Markdown(f'**{name}** *({axis.upper()} axis, Fundamental = {base_freq:.4f} Hz)*'))
        
        freqs, Pxx = spectogram(ts[axis])
        spectrograms.append((name, freqs, Pxx))
        plt.show()          # plt.savefig(f'x_axis_fft_{NFFT}.png')
    
    return spectrograms


def show_spectrogram_detail(spectrograms: list, axis: str, t: float):
    for name, freqs, Pxx in spectrograms:
        display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))
        i_window = spectrum_slice(freqs, Pxx, t)
        plt.show()           #plt.savefig(f'x_axis_fft_{NFFT}_at_{T_SEC}s.png')


def show_mms_peaks(spectrograms: list, axis: str, t: float):
    for name, freqs, Pxx in spectrograms:
        display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))
    
        i_window = window_idx(t)
        peaks = discovery.mms_peak_finder(Pxx.T[i_window])
        
        fig, ax = plt.subplots(1, 1, figsize=(15, 3))
        ax.grid(True)
        ax.plot(freqs, Pxx.T[i_window])
        ax.scatter(freqs[peaks], Pxx.T[i_window][peaks], marker='^', color='red')
        ax.set_xlabel('Frequency [Hz]')
        
        plt.show()


def show_harmonic_series(spectrograms: list, axis: str, t: float):
    # https://stackoverflow.com/questions/1982770/changing-the-color-of-an-axis
    for name, freqs, Pxx in spectrograms:
        display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))
    
        i_window = window_idx(t)
        h_series = discovery.harmonic_series_detection(freqs, Pxx.T[i_window], Fs, NFFT)
    
        # Find best (sum of harmonics' amplitudes in the largest)
        max_harmonic_amp_idx = np.argmax([
            sum([h[1] for h in s]) / len(s)
            for s in h_series
        ])
        best_harmonic_series = pd.DataFrame(
            h_series[max_harmonic_amp_idx],
            columns=['Frequency [Hz]', 'Amplitude [m/s^2]']
        )
        best_harmonic_series.index += 1
        display(tabulate(best_harmonic_series, headers='keys', tablefmt='html'))
    
        # Plot found harmonic series
        fig, ax = plt.subplots(1, 8, figsize=(30, 4))
        for i in range(8):
            s = h_series[i+1]
            if i == max_harmonic_amp_idx:
                ax[i].xaxis.label.set_color('red')
    
            ax[i].plot(freqs, Pxx.T[i_window])
            ax[i].scatter([x[0] for x in s], [x[1] for x in s], marker='^', color='red')
            ax[i].set_xlabel('Frequency [Hz]')
    
        plt.show()

def show_spectra_largest_amplitudes(spectrograms: list, axis: str, t: float):
    for name, freqs, Pxx in spectrograms:
        display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))

        i_window = window_idx(t)
        x_fundamental = get_max_frequency(freqs, Pxx, i_window)
        peaks = get_peaks(freqs, Pxx, i_window)
        
        display(Markdown(f'- *Fundamental frequency:* {x_fundamental} Hz'))
        display(tabulate(peaks.head(5), headers='keys', tablefmt='html'))


def compare_limited_specrograms(spectrograms: list, axis: str, t: float):
    fig, ax = plt.subplots(len(DATASET), 1, figsize=(20, 20), sharey=True)
    i = 0
    for name, ts in DATASET:
        signal = ts[axis].loc[t:t+NFFT/Fs].to_numpy()
        pxx = np.abs(np.fft.rfft(signal) / len(signal)) 
        freqs = np.fft.fftfreq(len(signal), d=1/Fs)[:len(pxx)]
        #ilast = len(freqs[freqs < F_LIMIT])
        
        ax[i].plot(freqs, pxx)
        ax[i].grid(True)
        ax[i].set_xlabel('Frequency [Hz]')
        ax[i].set_ylabel('Amplitude [m/s^2]')
        ax[i].set_xlim(0, F_LIMIT)
        ax[i].set_ylim(0, 0.4)
        ax[i].set_title(name)
        i += 1


def spectrogram_energy_left_cumulative(spectrograms: list, axis: str, t: float):
    fig, ax = plt.subplots(len(DATASET), 1, figsize=(20, 20), sharey=True)
    i = 0
    for name, ts in DATASET:
        signal = ts[axis].loc[t:t+NFFT/Fs].to_numpy()
        pxx = np.abs(np.fft.rfft(signal) / len(signal)) 
        freqs = np.fft.fftfreq(len(signal), d=1/Fs)[:len(pxx)]
        
        ax[i].plot(freqs, np.cumsum(pxx) / np.sum(pxx))
        ax[i].grid(True)
        ax[i].set_xlabel('Frequency [Hz]')
        ax[i].set_ylabel('Cumulative energy [%]')
        #ax[i].set_xlim(0, 10000)
        ax[i].set_title(name)
        i += 1

Compare mafaulda faults

In [ ]:
compare_limited_specrograms(DATASET, 'x', T_SEC)
plt.tight_layout()
plt.show()
No description has been provided for this image
In [ ]:
compare_limited_specrograms(DATASET, 'y', T_SEC)
plt.tight_layout()
plt.show()
No description has been provided for this image
In [ ]:
compare_limited_specrograms(DATASET, 'z', T_SEC)
plt.tight_layout()
plt.show()
No description has been provided for this image

Compare cumulative sums

In [ ]:
spectrogram_energy_left_cumulative(DATASET, 'x', T_SEC)
plt.tight_layout()
plt.show()
No description has been provided for this image
In [ ]:
spectrogram_energy_left_cumulative(DATASET, 'y', T_SEC)
plt.tight_layout()
plt.show()
No description has been provided for this image
In [ ]:
spectrogram_energy_left_cumulative(DATASET, 'z', T_SEC)
plt.tight_layout()
plt.show()
No description has been provided for this image

Spectrogram in X axis

In [ ]:
x_spectra = get_spectrograms(DATASET, 'x')

1_still.tsv (X axis, Fundamental = 18.4555 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

2_still.tsv (X axis, Fundamental = 20.3534 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

3_still.tsv (X axis, Fundamental = 21.9895 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

1_up.tsv (X axis, Fundamental = 18.5209 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

2_up.tsv (X axis, Fundamental = 20.4188 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

3_up.tsv (X axis, Fundamental = 22.0550 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

1_down.tsv (X axis, Fundamental = 17.2120 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

2_down.tsv (X axis, Fundamental = 20.3534 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

3_down.tsv (X axis, Fundamental = 21.3351 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

Spectrogram detail in X axis

In [ ]:
show_spectrogram_detail(x_spectra, 'x', T_SEC)

1_still.tsv (X axis @ 5s)

No description has been provided for this image

2_still.tsv (X axis @ 5s)

No description has been provided for this image

3_still.tsv (X axis @ 5s)

No description has been provided for this image

1_up.tsv (X axis @ 5s)

No description has been provided for this image

2_up.tsv (X axis @ 5s)

No description has been provided for this image

3_up.tsv (X axis @ 5s)

No description has been provided for this image

1_down.tsv (X axis @ 5s)

No description has been provided for this image

2_down.tsv (X axis @ 5s)

No description has been provided for this image

3_down.tsv (X axis @ 5s)

No description has been provided for this image

Peaks in frequency spectrum in X axis

  • MMS peak finder algorithm
In [ ]:
show_mms_peaks(x_spectra, 'x', T_SEC)

1_still.tsv (X axis @ 5s)

No description has been provided for this image

2_still.tsv (X axis @ 5s)

No description has been provided for this image

3_still.tsv (X axis @ 5s)

No description has been provided for this image

1_up.tsv (X axis @ 5s)

No description has been provided for this image

2_up.tsv (X axis @ 5s)

No description has been provided for this image

3_up.tsv (X axis @ 5s)

No description has been provided for this image

1_down.tsv (X axis @ 5s)

No description has been provided for this image

2_down.tsv (X axis @ 5s)

No description has been provided for this image

3_down.tsv (X axis @ 5s)

No description has been provided for this image

Harmonic series detection in X axis

In [ ]:
# show_harmonic_series(x_spectra, 'x', T_SEC)
In [ ]:
show_spectra_largest_amplitudes(x_spectra, 'x', T_SEC)

1_still.tsv (X axis @ 5s)

  • Fundamental frequency: 19.53125 Hz
f y 1x
0 19.53120.546024 1
1195.312 0.0096457610
2134.277 0.00462071 6.875
3214.844 0.0027503411
4412.598 0.0024323121.125

2_still.tsv (X axis @ 5s)

  • Fundamental frequency: 19.53125 Hz
f y 1x
0 19.53121.60475 1
1 51.26950.0283099 2.625
2195.312 0.010183 10
3273.438 0.0056699414
4332.031 0.0047059617

3_still.tsv (X axis @ 5s)

  • Fundamental frequency: 21.97265625 Hz
f y 1x
0 21.9727 1.49894 1
1314.941 0.011506 14.3333
2329.59 0.0085340815
3 97.6562 0.00731069 4.44444
4 9.765620.00419778 0.444444

1_up.tsv (X axis @ 5s)

  • Fundamental frequency: 17.08984375 Hz
f y 1x
0 17.08980.468103 1
1 70.80080.0112762 4.14286
2295.41 0.0067000617.2857
3144.043 0.00419564 8.42857
4 87.89060.00272608 5.14286

2_up.tsv (X axis @ 5s)

  • Fundamental frequency: 19.53125 Hz
f y 1x
0 19.53122.06101 1
1 78.125 0.0223562 4
2178.223 0.00997867 9.125
3307.617 0.0061318915.75
4314.941 0.0048611816.125

3_up.tsv (X axis @ 5s)

  • Fundamental frequency: 21.97265625 Hz
f y 1x
0 21.97271.00514 1
1195.312 0.0102214 8.88889
2292.969 0.0089147413.3333
3 12.207 0.00589523 0.555556
4163.574 0.00447052 7.44444

1_down.tsv (X axis @ 5s)

  • Fundamental frequency: 17.08984375 Hz
f y 1x
0 17.08980.0147896 1
1 53.71090.00494662 3.14286
2205.078 0.0019790812
3 83.00780.00169884 4.85714
4173.34 0.0015201310.1429

2_down.tsv (X axis @ 5s)

  • Fundamental frequency: 19.53125 Hz
f y 1x
0 19.53120.0323693 1
1 48.82810.00282413 2.5
2163.574 0.00189572 8.375
3361.328 0.0015244718.5
4 80.56640.00146388 4.125

3_down.tsv (X axis @ 5s)

  • Fundamental frequency: 21.97265625 Hz
f y 1x
0 21.97270.0823231 1
1 65.918 0.006334 3
2190.43 0.00178213 8.66667
3339.355 0.0014655515.4444
4283.203 0.0012035112.8889

Spectrogram in Y axis

In [ ]:
y_spectra = get_spectrograms(DATASET, 'y')

1_still.tsv (Y axis, Fundamental = 18.4555 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

2_still.tsv (Y axis, Fundamental = 20.3534 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

3_still.tsv (Y axis, Fundamental = 21.9895 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

1_up.tsv (Y axis, Fundamental = 18.5209 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

2_up.tsv (Y axis, Fundamental = 20.4188 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

3_up.tsv (Y axis, Fundamental = 22.0550 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

1_down.tsv (Y axis, Fundamental = 17.2120 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

2_down.tsv (Y axis, Fundamental = 19.5681 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

3_down.tsv (Y axis, Fundamental = 21.3351 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

Spectrogram detail in Y axis

In [ ]:
show_spectrogram_detail(y_spectra, 'y', T_SEC)

1_still.tsv (Y axis @ 5s)

No description has been provided for this image

2_still.tsv (Y axis @ 5s)

No description has been provided for this image

3_still.tsv (Y axis @ 5s)

No description has been provided for this image

1_up.tsv (Y axis @ 5s)

No description has been provided for this image

2_up.tsv (Y axis @ 5s)

No description has been provided for this image

3_up.tsv (Y axis @ 5s)

No description has been provided for this image

1_down.tsv (Y axis @ 5s)

No description has been provided for this image

2_down.tsv (Y axis @ 5s)

No description has been provided for this image

3_down.tsv (Y axis @ 5s)

No description has been provided for this image

Peaks in frequency spectrum in Y axis

In [ ]:
show_mms_peaks(y_spectra, 'y', T_SEC)

1_still.tsv (Y axis @ 5s)

No description has been provided for this image

2_still.tsv (Y axis @ 5s)

No description has been provided for this image

3_still.tsv (Y axis @ 5s)

No description has been provided for this image

1_up.tsv (Y axis @ 5s)

No description has been provided for this image

2_up.tsv (Y axis @ 5s)

No description has been provided for this image

3_up.tsv (Y axis @ 5s)

No description has been provided for this image

1_down.tsv (Y axis @ 5s)

No description has been provided for this image

2_down.tsv (Y axis @ 5s)

No description has been provided for this image

3_down.tsv (Y axis @ 5s)

No description has been provided for this image

Harmonic series detection in Y axis

In [ ]:
# show_harmonic_series(y_spectra, 'y', T_SEC)
In [ ]:
show_spectra_largest_amplitudes(y_spectra, 'y', T_SEC)

1_still.tsv (Y axis @ 5s)

  • Fundamental frequency: 19.53125 Hz
f y 1x
0 19.53120.823792 1
1117.188 0.0125524 6
2 73.24220.00764154 3.75
3214.844 0.0057619311
4166.016 0.00391798 8.5

2_still.tsv (Y axis @ 5s)

  • Fundamental frequency: 19.53125 Hz
f y 1x
0 19.53121.63041 1
1 53.71090.0175563 2.75
2183.105 0.0129831 9.375
3292.969 0.0072152615
4266.113 0.0056260913.625

3_still.tsv (Y axis @ 5s)

  • Fundamental frequency: 21.97265625 Hz
f y 1x
0 21.97271.86651 1
1109.863 0.020548 5
2393.066 0.0126337 17.8889
3361.328 0.009655 16.4444
4131.836 0.00747923 6

1_up.tsv (Y axis @ 5s)

  • Fundamental frequency: 17.08984375 Hz
f y 1x
0 17.08980.684795 1
1104.98 0.00932373 6.14286
2 75.68360.00481567 4.42857
3275.879 0.0030466216.1429
4305.176 0.0022769817.8571

2_up.tsv (Y axis @ 5s)

  • Fundamental frequency: 19.53125 Hz
f y 1x
0 19.53121.96031 1
1178.223 0.0153813 9.125
2185.547 0.00966599 9.5
3158.691 0.00717196 8.125
4244.141 0.0063370312.5

3_up.tsv (Y axis @ 5s)

  • Fundamental frequency: 21.97265625 Hz
f y 1x
0 21.97271.66467 1
1197.754 0.017065 9
2109.863 0.0124474 5
3 75.68360.007274673.44444
4 97.65620.006595814.44444

1_down.tsv (Y axis @ 5s)

  • Fundamental frequency: 2.44140625 Hz
f y 1x
0 2.441410.0271932 1
1114.746 0.00266754 47
2476.074 0.0018169 195
3173.34 0.00161734 71
4522.461 0.00147466 214

2_down.tsv (Y axis @ 5s)

  • Fundamental frequency: 19.53125 Hz
f y 1x
0 19.53120.0651357 1
1 51.26950.00371688 2.625
2131.836 0.00216978 6.75
3102.539 0.00176165 5.25
4549.316 0.0015881128.125

3_down.tsv (Y axis @ 5s)

  • Fundamental frequency: 21.97265625 Hz
f y 1x
0 21.97270.0480764 1
1 85.44920.00466419 3.88889
2244.141 0.0020465 11.1111
3300.293 0.0017056 13.6667
4280.762 0.0015513712.7778

Spectrogram in Z axis

In [ ]:
z_spectra = get_spectrograms(DATASET, 'z')

1_still.tsv (Z axis, Fundamental = 13.4817 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

2_still.tsv (Z axis, Fundamental = 20.3534 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

3_still.tsv (Z axis, Fundamental = 21.9895 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

1_up.tsv (Z axis, Fundamental = 0.0654 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

2_up.tsv (Z axis, Fundamental = 0.0654 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

3_up.tsv (Z axis, Fundamental = 22.0550 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

1_down.tsv (Z axis, Fundamental = 0.0654 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

2_down.tsv (Z axis, Fundamental = 0.0654 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

3_down.tsv (Z axis, Fundamental = 0.0654 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz
No description has been provided for this image

Spectrogram detail in Z axis

In [ ]:
show_spectrogram_detail(z_spectra, 'z', T_SEC)

1_still.tsv (Z axis @ 5s)

No description has been provided for this image

2_still.tsv (Z axis @ 5s)

No description has been provided for this image

3_still.tsv (Z axis @ 5s)

No description has been provided for this image

1_up.tsv (Z axis @ 5s)

No description has been provided for this image

2_up.tsv (Z axis @ 5s)

No description has been provided for this image

3_up.tsv (Z axis @ 5s)

No description has been provided for this image

1_down.tsv (Z axis @ 5s)

No description has been provided for this image

2_down.tsv (Z axis @ 5s)

No description has been provided for this image

3_down.tsv (Z axis @ 5s)

No description has been provided for this image

Peaks in frequency spectrum in Z axis

In [ ]:
show_mms_peaks(z_spectra, 'z', T_SEC)

1_still.tsv (Z axis @ 5s)

No description has been provided for this image

2_still.tsv (Z axis @ 5s)

No description has been provided for this image

3_still.tsv (Z axis @ 5s)

No description has been provided for this image

1_up.tsv (Z axis @ 5s)

No description has been provided for this image

2_up.tsv (Z axis @ 5s)

No description has been provided for this image

3_up.tsv (Z axis @ 5s)

No description has been provided for this image

1_down.tsv (Z axis @ 5s)

No description has been provided for this image

2_down.tsv (Z axis @ 5s)

No description has been provided for this image

3_down.tsv (Z axis @ 5s)

No description has been provided for this image

Harmonic series detection in Z axis

In [ ]:
# show_harmonic_series(z_spectra, 'z', T_SEC)
In [ ]:
show_spectra_largest_amplitudes(z_spectra, 'z', T_SEC)

1_still.tsv (Z axis @ 5s)

  • Fundamental frequency: 14.6484375 Hz
f y 1x
0 14.64840.0702139 1
1 24.41410.00735616 1.66667
2104.98 0.00376746 7.16667
3163.574 0.0026874611.1667
4231.934 0.0020599 15.8333

2_still.tsv (Z axis @ 5s)

  • Fundamental frequency: 19.53125 Hz
f y 1x
0 19.53120.158409 1
1217.285 0.045319 11.125
2 31.73830.0148314 1.625
3141.602 0.00963871 7.25
4334.473 0.0046994517.125

3_still.tsv (Z axis @ 5s)

  • Fundamental frequency: 21.97265625 Hz
f y 1x
0 21.97270.275604 1
1 65.918 0.00763671 3
2292.969 0.0047130813.3333
3205.078 0.00322265 9.33333
4183.105 0.00215798 8.33333

1_up.tsv (Z axis @ 5s)

  • Fundamental frequency: 14.6484375 Hz
f y 1x
0 14.64840.0538022 1
1 26.85550.00976307 1.83333
2 61.03520.00524637 4.16667
3236.816 0.0040075616.1667
4393.066 0.0020804126.8333

2_up.tsv (Z axis @ 5s)

  • Fundamental frequency: 19.53125 Hz
f y 1x
0 19.53120.296585 1
1117.188 0.0501181 6
2178.223 0.0220361 9.125
3107.422 0.0153412 5.5
4205.078 0.008809510.5

3_up.tsv (Z axis @ 5s)

  • Fundamental frequency: 21.97265625 Hz
f y 1x
0 21.97270.211261 1
1175.781 0.006943728
2148.926 0.003578596.77778
3 73.24220.002992983.33333
4107.422 0.002058584.88889

1_down.tsv (Z axis @ 5s)

  • Fundamental frequency: 2.44140625 Hz
f y 1x
0 2.441410.0193371 1
1 92.7734 0.0044802 38
2214.844 0.00232566 88
3173.34 0.00186004 71
4371.094 0.00147204 152

2_down.tsv (Z axis @ 5s)

  • Fundamental frequency: 2.44140625 Hz
f y 1x
0 2.441410.0535651 1
1 65.918 0.00347521 27
2144.043 0.0028824 59
3153.809 0.00226953 63
4187.988 0.00185768 77

3_down.tsv (Z axis @ 5s)

  • Fundamental frequency: 14.6484375 Hz
f y 1x
0 14.64840.0340945 1
1 63.47660.00477927 4.33333
2256.348 0.0028789317.5
3180.664 0.0019068412.3333
4 90.332 0.00168929 6.16667

Histogram

In [ ]:
axis = ['x', 'y', 'z']
for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    ts[axis].hist(figsize=(10, 5), grid=True, bins=50)
    plt.show()

1_still.tsv

No description has been provided for this image

2_still.tsv

No description has been provided for this image

3_still.tsv

No description has been provided for this image

1_up.tsv

No description has been provided for this image

2_up.tsv

No description has been provided for this image

3_up.tsv

No description has been provided for this image

1_down.tsv

No description has been provided for this image

2_down.tsv

No description has been provided for this image

3_down.tsv

No description has been provided for this image
In [ ]:
axis = ['x', 'y', 'z']
for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    ts[axis].boxplot(figsize=(10, 5))
    plt.show()

1_still.tsv

No description has been provided for this image

2_still.tsv

No description has been provided for this image

3_still.tsv

No description has been provided for this image

1_up.tsv

No description has been provided for this image

2_up.tsv

No description has been provided for this image

3_up.tsv

No description has been provided for this image

1_down.tsv

No description has been provided for this image

2_down.tsv

No description has been provided for this image

3_down.tsv

No description has been provided for this image

Orbitals of all cross sections

In [ ]:
for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    fig, ax = plt.subplots(1, 3, figsize=(20, 4))

    for i, col in enumerate([('x', 'y'), ('x', 'z'), ('y', 'z')]):
        ax[i].scatter(ts[col[0]], ts[col[1]], s=1)
        ax[i].grid(True)
        ax[i].set_xlabel(col[0].upper())
        ax[i].set_ylabel(col[1].upper())
        ax[i].grid(True)
    plt.show()       # plt.savefig('orbitals.png')

1_still.tsv

No description has been provided for this image

2_still.tsv

No description has been provided for this image

3_still.tsv

No description has been provided for this image

1_up.tsv

No description has been provided for this image

2_up.tsv

No description has been provided for this image

3_up.tsv

No description has been provided for this image

1_down.tsv

No description has been provided for this image

2_down.tsv

No description has been provided for this image

3_down.tsv

No description has been provided for this image

Orbitals of 1x harmonic frequency

In [ ]:
x_spectra_by_name = {spec[0]: spec for spec in x_spectra}
y_spectra_by_name = {spec[0]: spec for spec in y_spectra}
z_spectra_by_name = {spec[0]: spec for spec in z_spectra}
t = 5
space = 5

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    fig, ax = plt.subplots(1, 3, figsize=(20, 4))

    name, freqs, Pxx = x_spectra_by_name[name]
    x_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))
    name, freqs, Pxx = y_spectra_by_name[name]
    y_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))

    name, freqs, Pxx = z_spectra_by_name[name]
    z_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))

    try:
        ts['x_1x'] = butter_bandpass_filter(ts['x'], x_fundamental - space, x_fundamental + space, Fs)
        ts['y_1x'] = butter_bandpass_filter(ts['y'], y_fundamental - space, y_fundamental + space, Fs)
        ts['z_1x'] = butter_bandpass_filter(ts['z'], z_fundamental - space, z_fundamental + space, Fs)
    except ValueError:
        continue
    
    for i, col in enumerate([('x_1x', 'y_1x'), ('x_1x', 'z_1x'), ('y_1x', 'z_1x')]):
        ax[i].scatter(ts[col[0]], ts[col[1]], s=1)
        ax[i].grid(True)
        ax[i].set_xlabel(col[0].upper())
        ax[i].set_ylabel(col[1].upper())
        ax[i].grid(True)
    
    plt.show()       # plt.savefig('orbitals_1x.png')

1_still.tsv

No description has been provided for this image

2_still.tsv

No description has been provided for this image

3_still.tsv

No description has been provided for this image

1_up.tsv

No description has been provided for this image

2_up.tsv

No description has been provided for this image

3_up.tsv

No description has been provided for this image

1_down.tsv

2_down.tsv

3_down.tsv

No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
In [ ]:
t = 5
space = 8

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))

    name, freqs, Pxx = x_spectra_by_name[name]
    x_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))
    name, freqs, Pxx = y_spectra_by_name[name]
    y_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))

    name, freqs, Pxx = z_spectra_by_name[name]
    z_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))

    try:
        x = butter_bandpass_filter(ts['x'], x_fundamental - space, x_fundamental + space, Fs)
        y = butter_bandpass_filter(ts['y'], y_fundamental - space, y_fundamental + space, Fs)
        z = butter_bandpass_filter(ts['z'], z_fundamental - space, z_fundamental + space, Fs)
    except ValueError:
        continue

    ax = plt.figure().add_subplot(projection='3d')
    ax.scatter(x, y, z, zdir='z', s=1, color='navy')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')
    ax.set_xlim(-2, 2)
    ax.set_ylim(-2, 2)
    ax.set_zlim(-2, 2)
    ax.zaxis.labelpad = -0.7
    plt.show()

1_still.tsv

No description has been provided for this image

2_still.tsv

No description has been provided for this image

3_still.tsv

No description has been provided for this image

1_up.tsv

No description has been provided for this image

2_up.tsv

No description has been provided for this image

3_up.tsv

No description has been provided for this image

1_down.tsv

2_down.tsv

3_down.tsv

No description has been provided for this image